ParamSe=zeros(NZ,1); %Standard error of estimated coefficients of Z variables
MeanSE=zeros(NV,1); %Standard error of Mean
StdSE=zeros(NV,1);  %  " of Std Devs
CorrSE=zeros(NV,NV);  % " of Correlation matrix
FreqSE=zeros(NV,NBins); % " of Freq in each bin
MidSE=zeros(NV,NBins); % " of Midpoint for each bin, (std err should be zero)

paramhold=zeros(NZ,NReps);
llhold=zeros(1,NReps);
mnhold=zeros(NV,NReps);
stdhold=zeros(NV,NReps);
corrhold=zeros(NV-2,NV-2,NReps);
if discountfun=="QuasiHyper"
    corrhold = zeros(NV-3, NV-3, NReps);
end
freqhold=zeros(NV,NBins,NReps);
midhold=zeros(NV,NBins,NReps);

XMAT_Original=XMAT;
PID_Original=XMAT(:,1);
CSID_Original=XMAT(:,2);
clear i r

llorig = finalLL;
qualitytime0 = qualitytime;
wlifetime0 = wlifetime;
infratime0 = infratime;
jobstime0 = jobstime;
buffertime0 = buffertime;
swtime0 = swtime;
captime0 = captime;


for xi=1:NReps
    FirstRun=0;
    disp('Bootstrap Sample #');
    disp(xi);
    disp(datetime);
    disp('');
    btsample=randi(NP,[NP,1]);
    XMAT=[];
    qualitytime = [];
    wlifetime = [];
    infratime = [];
    jobstime = [];
    buffertime = [];
    swtime= []; 
    captime= [];
    for r=1:NP;
       thisp=btsample(r,1);
       rowsr = PID_Original==thisp;
       thisx=XMAT_Original(rowsr,:);
       thisx(:,1)=r.*ones(size(thisx,1),1); %creates the ID number sequentially
       XMAT=[XMAT ; thisx];   
       captime = [captime; captime0(rowsr,:)];
       swtime = [swtime; swtime0(rowsr,:)];
       buffertime = [buffertime; buffertime0(rowsr,:)];
       qualitytime = [qualitytime; qualitytime0(rowsr,:)];
       wlifetime = [wlifetime; wlifetime0(rowsr,:)];
       infratime = [infratime; infratime0(rowsr,:)];
       jobstime = [jobstime; jobstime0(rowsr,:)];
    end
    XMAT(:,2)=CSID_Original;
    if discountfun=="QuasiHyper"
        CreateDrawsDisWtpProbsQH;
        CreateZQH;
    else
        CreateDrawsDisWtpProbs;
        CreateZ;
    end
    if YesGPU==1       
        PROBS=gpuArray(PROBS);
        BETAS=gpuArray(BETAS);
        Z=gpuArray(Z);
    end
    param=paramhat_orig; %StartB; %Trying to start from estimate of real dataset
    options=optimset('LargeScale','off','Display','final','GradObj','on',...
    'MaxFunEvals',10000,'MaxIter',MAXITERS,'TolX',PARAMTOL,'TolFun',LLTOL,'DerivativeCheck','off');
    [paramhat,fval,exitflag]=fminunc(@flexll,param,options);
    if YesGPU==1;
        PROBS=gather(PROBS);
        BETAS=gather(BETAS);
        Z=gather(Z);
    end;
    paramhold(:,xi)=paramhat;
    llhold(1,xi) = finalLL;
    [mnhold(:,xi),stdhold(:,xi),cc,freqhold(:,:,xi),midhold(:,:,xi)]=stats(paramhat,NBins);
    if discountfun=="QuasiHyper"
        corrhold(:,:,xi)=corrcov(cc(4:end,4:end));
    else
        corrhold(:,:,xi)=corrcov(cc(3:end,3:end));
    end
    %clear f;
    save(filename);
end
finalLL = llorig;
qualitytime = qualitytime0;
wlifetime = wlifetime0;
infratime = infratime0;
jobstime = jobstime0;
buffertime = buffertime0;
swtime = swtime0;
captime = captime0;
XMAT=XMAT_Original;
clear llorig XMAT_Original qualitytime0 wlifetime0 infratime0 jobstime0 buffertime0 swtime0 captime0;
ParamSE=std(paramhold,0,2);
MeanSE=std(mnhold,0,2);
StdSE=std(stdhold,0,2);
CorrSE=std(corrhold,0,3); 
FreqSE=std(freqhold,0,3); 
MidSE=std(midhold,0,3);

disp('Final Log-likelihood');
disp(' ');
disp(finalLL);
disp(' ');
disp('Utility Parameters');
disp(' ');
disp('                    Mean                   Std Dev');
disp('              ------------------   -----------------------');
disp('                 Est     SE            Est         SE');

for r=1:NV
        fprintf('%-10s %10.4f %10.4f %10.4f %10.4f\n', NAMES{r,1}, [MeanEst(r,1) MeanSE(r,1) StdEst(r,1) StdSE(r,1)]);
end
if discountfun=="QuasiHyper"
    disp(' ');
    disp('Correlations of WTPs');
    disp(corrcov(CovMatEst(4:end,4:end)));
    disp('T-stats on Correlations');
    disp('Diagonal has Inf because diagonals of correlation matrix are always 1.');
    disp(corrcov(CovMatEst(4:end,4:end))./CorrSE);
    disp(' ');
else
    disp(' ');
    disp('Correlations of WTPs');
    disp(corrcov(CovMatEst(3:end,3:end)));
    disp('T-stats on Correlations');
    disp('Diagonal has Inf because diagonals of correlation matrix are always 1.');
    disp(corrcov(CovMatEst(3:end,3:end))./CorrSE);
    disp(' ');
end
Lm=zeros(NV,1);
Um=zeros(NV,1);
Ls=zeros(NV,1);
Us=zeros(NV,1);
Lp=zeros(NZ,1);
Up=zeros(NZ,1);
for(ii=1:(NZ))
    Lp(ii,1)=quantile(paramhold(ii,:),0.025);
    Up(ii,1)=quantile(paramhold(ii,:),0.975);
end
for(ii=1:NV)
    Lm(ii,1)=quantile(mnhold(ii,:), 0.025);
    Um(ii,1)=quantile(mnhold(ii,:), 0.975);
    Ls(ii,1)=quantile(stdhold(ii,:), 0.025);
    Us(ii,1)=quantile(stdhold(ii,:), 0.975);
end

BootCIparam = horzcat(paramhat_orig, ParamSE, Lp,Up);
tab = dataset({BootCIparam 'Estimate', 'SE', 'Lower', 'Upper'},...
    'obsnames', {});
tab = dataset2table(tab);
nm = horzcat('ParamResults.',modelname, '.csv');
writetable(tab, nm);

%disp('First Row of ParamResults.csv contains Price coefficent if in preference space')

BootCImean = horzcat(MeanEst, MeanSE, Lm,Um, StdEst, StdSE, Ls, Us);

Results = dataset({BootCImean 'MeanEst', 'MeanSE', 'LowerMean', 'UpperMean', 'SDEst', 'SDSE', 'LowerSE', 'UpperSE'},...
    'obsnames', NAMES);

Results = dataset2table(Results);
nm = horzcat('DistrResults.',modelname, '.csv');
writetable(Results, nm,'WriteRowNames',true);

disp(Results);
save(filename)

